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EXECUTIVE  SUMMARY 


This  report  discusses  the  integration  of  a  novel  finite  element  for  analyzing  composite 
material  structures  into  a  commercially  available  finite  element  software  package.  This  element  is 
a  constitutive  equation  based  formulation  that  predicts  individual  layer  strains  and  stresses  in  a 
computationally  efficient  manner.  Implementing  this  element  into  a  commercial  code  will  improve 
its  efficiency  for  predicting  the  structural  response  of  composite  material  structures. 

The  work  discussed  in  this  report  involved  utilizing  algorithms  previously  developed  at 
Southwest  Research  Institute  (SwRI)  that  are  based  on  the  constitutive  equations  and  integrating 
them  with  the  ANSYS  finite  element  code.  An  example  problem  using  this  method  is  presented  that 
predicts  the  laminate  free  edge  stress  distribution  that  is  important  for  examining  potential 
deiaminations  at  design  features  such  as  holes,  cutouts,  and  joints. 

1  he  results  of  this  study  show  that  the  element  is  accurate  for  determining  the  stress 
distribution  at  the  laminate  tree  edge.  Recommendations  for  enhancements  to  the  element  are  made 
to  improve  its  modeling  efficiency  even  further. 
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1.0  lN'TROinXTlON 


This  report  discusses  \iic  integration  of  a  novel  finite  element  tor  analyzing  composite' 
materia!  structures  into  a  commercioi  finite  element  software  code.  The  work  described  within  was 
Task  3.0  of  the  L:.S.  Govet.unent  contract  N00 1 67-95-C -0064.  and  is  intended  as  a  stand-alone 
document.  The  government  report  numbered  CD/NSWC  MSSPO  (102)  CR- 95/02  and  entitled 
"Development  of  Large.  Thick  Laminate  Structures  for  Offshore  Applications'",  discusses  the  rest 
of  the  tasks  associated  with  this  program. 

The  thrust  of  this  effort  was  to  develop  a  computationally  effective  method  to  anah/.c  the 
through-thickness  stress  and  strain  distribution  of  thick  composite  parts  as  applied  to  offshore 
structures.  However  this  method  can  used  on  a  variety  of  composite  material  applications  for 
structures. 

This  report,  is  organized  into  the  following  sections:  background  information  and  the  rationale 
for  using  this  method  are  presented  in  Section  2.0.  Section  3.0  discusses  the  technical  approach  for 
integrating  the  constitutive  based  model  into  a  con  inetcial  finite  element  program.  The  results  of 
an  example  problem  based  on  the  free  edge  stress  distribution  problem  are  presented  in  Section  4.0. 
Section  5.0  discusses  the  conclusions  of  using  this  method.  Recommendations  for  future  efforts  are 
discussed  in  Section  6.0.  Appendix  A  presents  the  source  code  listing  for  integrating  the  model  into 
the  commercial  finite  element  code.  Results  from  the  example  problem  are  presented  in 
Appendix  B. 


2.0  B.U  K(.ROIM) 


The  application  of  composite  structures  tor  otYshore  structures  will  require  parts  that  arc  thick 
(relative  to  typical  aerospace  structures),  and  ha\e  loads  that  are  applied  through  the  thickness  plane 
of  the  structure  Particular  attention  must  be  paid  to  potential  delamination  of  the  composite  material 
near  design  details  such  as  holes,  cutouts,  and  both  compositc-to-coinpositc  and  steel-to-compositc 
joint  regions. 

l'herefore  an  analytical  method  that  can  determine  the  structural  response  of  these  laminated 
parts  must  be  used  to  ensure  their  soundness  when  placed  into  service.  Due  to  the  complex  geo  metre 
and  loading  conditions  acting  on  this  parts,  the  finite  element  method  must  be  employed  to 
adequately  analyze  their  structural  behavior.  Moreover,  the  method  used  mast  be  computationally 
efficient  to  permit  analyzing  these  structures  in  a  timely  and  cost-effective  manner. 

There  are  many  methods  available  to  describe  the  mechanical  response  of  layered  composite 
materials  Typically  these  methods  are  based  on  classical  laminated  plate  theory  with  estimates  of 
the  bulk  mechanical  response  for  thick  composite  sections.  Classical  laminate  plate  theory  is  v  ery 
effective  in  modeling  the  response  of  thin  laminates  when  there  are  no  significant  loads  applied  in 
the  thickness  direction  of  the  laminate.  Howev  er,  when  this  is  not  the  case,  some  other  material 
description  must  be  used  to  determine  the  through-thickness  mechanical  response. 

One  common  approach  to  describe  the  through  thickness  material  properties  is  to  replace  the 
layer  by  layer  characteristics  with  an  average  set  of  orthotropic  constitutiv  e  properties  that  are  used 
to  model  many  layers  of  composite.  This  method  is  adequate  for  large  sections  of  composite,  but 
fail  to  determine  the  layer-by-layer  behavior  that  is  critical  information  near  design  details.  One 
means  to  get  around  this  problem  is  to  employ  “Global-Local'*  modeling,  where  a  v  ery  refined  layer- 
by-layer  finite  model  mesh  is  used  around  the  local  design  details,  and  a  coarser  mesh  with  effective 
properties  is  used  for  the  bulk  response  of  the  structure.  Two  major  problems  with  this  method  are 
that  the  analyst  must  know  beforehand  what  design  features  are  important  for  the  bulk  mechanical 
response  of  the  structure,  and  if  there  are  many  design  features  ef  interest  then  the  refined  layer-by- 
layer  mesh  portions  of  the  model  can  grow  to  the  point  where  computation  becomes  prohibitively 
expensive.  It  is  clear  that  what  is  needed  is  a  reliable  three  dimensional  material  description  that  can 
be  implemented  into  a  finite  element  that  captures  the  essential  features  of  a  layered  composite 
without  the  computational  expense  of  a  layer-by-layer  approach. 

There  are  many  methods  proposed  in  the  literature  for  the  description  of  the  three 
dimensional  “effective"  properties  of  a  laminated  composite.  How'evcr,  these  methods  generally 
have  very  complex  mathematical  descriptions  of  these  properties  and  require  extensive  numerical 
matrix  manipulations.  These  characteristics  result  in  finite  element  solutions  that  are 
computationally  prohibitive  for  modeling  composite  structures  of  any  great  size  or  detail. 

A  method  that  overcomes  these  obstacles  is  a  generalized  averaging  procedure  for  the 
description  of  the  mechanical  properties  of  a  layered  composite  that  was  proposed  by 
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thristcnsen  [  1  |  l  lie  main  point  of  this  procedure  is  the  separation  of  The  liber-dominated  properties 
from  the  matrix  dominated  properties.  The  matrix  dominated  properties  are  described  In  two 
independent  constants  taat  are  the  result  of  this  averaging  procedure,  and  are  used  in  a  manner  that 
is  similar  to  Hooke’s  Taw.  The  fiber  dominated  properties  are  then  averaged  into  one  constant  that 
the  provides  an  additive  term  to  the  constitutive  formulation  in  tire  direction  of  the  fibers  in  that  layer 
of  the  composite. 

This  constitutive  equation  method  was  adapted  into  a  displacement-based  three  dimensional 
finite  element  code  by  Patton  [2].  This  code,  which  will  be  referred  io  as  the  SwRI  code,  is  a 
"university-grade"  code,  meaning  it  does  not  have  a  pro-  and  post  processor  and  has  limited 
functionality,  in  order  to  utilize  this  method  in  a  cost-effective  maimer,  it  is  necessary  to  implement 
it  into  a  commercial  grade  FT  A  code. 


3.0  TKCMMCAI.  APPROACH 


T  he  technical  approach  of  this  effort  was  to  use  portions  of  ihe  SwRI  code  specific  to  the 
single  element  formulation  ami  integrate  it  with  a  commeieial  finite  element  program  There  arc 
several  factors  that  must  be  considered  before  deciding  01  which  software  package  to  use. 

1  he  most  important  characteristic  for  choosing  which  finite  clement  code  to  integrate  the 
constitutive  equation  based  element  is  that  the  commercial  code  have  the  capability  to  accept  user 
written  elements.  Also  the  present  formulation  of  the  constitutive  equation  based  element  requires 
that  the  finite  element  program  have  a  displacement  based  wavefront  sober.  This  is  the  most 
common  means  of  solving  finite  element  models,  so  there  are  a  large  number  of  codes  a\  ailable. 
Another  requirement  is  that  the  program  preferably  have  an  integrated  pre-  and  post  processor,  or 
at  least  be  able  to  interface  with  a  third-party  pre-  and  post  processing  program  such  as  PAIR  AN 
or  IDT  AS.  A  highly  desiraole  requirement  is  that  the  program  is  already  in  place  within  the 
organization  and  has  an  experienced  user  base. 

There  are  a  number  of  finite  element  codes  that  meet  these  requiiemems.  NAS  1  RAN. 
ANSYS,  and  ABAQUS  among  the  better  known  products.  There  is  no  one  program  that  is  the  best 
choice,  the  decision  being  based  on  cost,  familiarity,  and  of  co  arse  whether  or  no;  if  one  alread\ 
possess  one  of  these  codes.  The  host  choice  for  the  Marine  Technology  Department  at  SwRI  was 
ANSYS  5.1  coded  to  run  a  HP  700  series  UNIX  Workstation.  It  should  he  noted  that  the  user 
element  interlace  capability  is  not  available  on  the  PC  -DCS  version  of  this  software  package. 

An  interface  to  the  ANSYS  code  is  provided  by  the  use  of  User  Programmable  Features 
(UPFs)  [3J.  UPFs  are  a  set  of  FORTRAN  routines  that  allow  the  user  to  write  custom  elements, 
loading  conditions,  material  properties,  and  other  user  specific  applications. 

This  particular  application  required  the  use  of  the  user  element  subroutines.  These 
subroutines  provide  the  interface  between  the  user  element  and  the  ANSYS  code.  The  user  is 
responsible  for  developing  the  entire  element  formulation  including  the  load  vector,  stiffness  matrix, 
and  material  matrix.  The  user  must  also  address  issues  with  ANSYS  variable  passing,  input/output 
switches,  and  any  options  for  the  element. 

The  portions  specific  to  single  element  formulations  in  the  SwRI  code  were  used  in  the 
ANSYS  user  routine.  Because  the  SwRI  code  was  developed  on  a  DOS  based  PC  with  Microsoft 
FORT  RAN  V5.  some  of  the  code  had  to  be  changed  to  successfully  port  it  to  the  HP  Workstation. 
The  source  code  listing  for  this  element  is  presented  in  Appendix  A. 
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An  example  problem  was  tested  to  determine  the  aeeuraey  and  utility  ol  implementing  the 
constitutive  formulation  with  ANSYS.  The  gi’ometiy.  material  properties  and  boundary  conditions 
for  the  example  problem  were  identical  to  the  one  solved  by  Pipes  and  I’agano  j4].  This  problem 
is  a  free  edge  stiess  problem  that  consists  of  a  4-ply  laminate  under  uniform  axial  extension.  Phis 
is  one  of  the  nu.st  demanding  solutions  for  laminate  analyses,  in  that  there  are  interlaminai  stresses 
at  the  free  edge  of  the  laminate  that  arc  not  accounted  for  by  classical  laminate  theory.  Additionally . 
this  same  problem  was  solved  by  the  existing  SwRI  code  and  compared  to  the  results  from  the 
ANSYS  model. 

4.1  Model  Description 

The  geometry  and  corresponding  finite  element  mesh  was  created  in  the  interactive  “mode" 
of  ANSYS.  which  is  a  graphical  interface.  This  model  is  a  llat  plate  that  is  1 .2  inches  long,  by  0.8 
inches  wide  by  0.2  inches  thick.  As  shown  in  f  igure  1.  it  is  a  +45/-45Z-  45  A  45  composite  structure, 
The  finite  element  model  consists  of  %  8-noded  brick  elements  and  is  shown  in  f  igure  2.  1  here  are 
two  symmetry  planes,  one  along  the  x-axis.  and  one  along  the  v-axis.  The  loading  is  a  uni-axial 
extension  applied  by  displacements  along  one  surface  of  the  plate  in  the  x-direc'ion.  The  value  of 
these  displacements  w  as  0.075  inch  at  each  node  of  the  surface. 

The  material  properties  for  this  problem  are  shown  in  Table  1 .  1  hese  values  are  input  using 
the  ANSYS  MP  command  via  the  command  line. 


TABLE  1-  MATERIAL  PROPERTIES 


j  Material  Property 

ANSYS 

Variable 

Value  | 

Young's  Modulus  along  liber  direction 

f.X  1 

2. 1  x  1()7  psi  | 

Young's  Modulus  perpendicular  to  fiber  direction 

f:Y 

Possions  Ratio 

NUXY 

021 

Shear  Modulus  in  xy  plane 

GXY 

0.85x10"  psi 

Shear  Modulus  in  xz  plane 

GXZ 

0.0 

Density 

RHO 

0.0 

Coefficient  of  thermal  expansion  in  x  direction 

ALPMAX 

0.0 

Coefficient  of  thermal  expansion  in  y  direction 

ALPI1AY 

0.0 

Coefficient  of  thermal  expansion  in  z  direction 

ALPHAZ 

0.0 

1  he  individual  ply  layer  properties  are  input  using  ANSYS  real  constants.  It  was  found  that 
the  parameters  for  the  layers  musT  be  entered  with  two  material  properties  through  the  lamina 


FIGURE  1  -  Ply  Orientation  of  Plate 


FIGURE  2  -  Finite  Element  Model 
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thickness.  This  requires  that  there  are  two  elements  through  the  thickness,  each  with  a  different 
ANSYS  material  property  ID  number.  However,  the  material  properties  are  identical  for  both 
material  sets.  It  is  suspected  that  this  condition  is  caused  by  the  formulation  of  the  material  matrix 
in  the  SwRI  code  and  will  eventually  need  to  be  rectified.  Therefore  a  separate  real  constant  is 
required  for  each  set  of  elements  through  the  laminate  thickness. 

The  real  constant  data  consists  of  the  number  of  layers  through  the  element,  layer  material 
numbers,  fiber  orientations,  and  layer  thickness.  This  data  is  entered  using  the  ANSYS  commands 
R  and  RMORE  via  the  command  line.  The  format  for  entering  the  real  constants  is  based  on  an 
existing  ANSYS  composite  element,  SOLID  46,  and  is  shown  below: 

R  ,Nl,NLR2,R3,R4,R5,R6 

Where:  N 1  is  the  real  constant  ID  set 

NL  is  the  number  of  layers  in  this  set 
R2  -  R6  are  not  used 

RMORE,;?  7,  R8,  R  9,  R 1 0,  R 1 1 ,  R 12 

Where:  R7  -  R1 2  are  not  used 

RMORE  JU3.R14.R15.R16.R17.R18 

Where:  R13  is  the  material  ID  for  the  first  layer  in  this  set 

R14  is  the  fiber  direction  from  element  x-axis  for  first  layer  in  this  set 
R1 5  is  the  lamina  thickness  of  the  first  layer  in  this  set 
R1 6  is  the  material  ID  for  the  second  layer  in  this  set 
R1 7  is  the  fiber  direction  from  element  x-axis  for  second  layer  in  this  set 
R18  is  the  lamina  thickness  of  the  second  layer  in  this  set 

The  command  line  data  entry  for  this  problem  is  as  follows: 

(First  real  constant  set,  layers  1-2) 

R  , 1  /  2 ,  ,  ,  ,  ,  , 

RMORE, , , , , , 

RMORE, 1,45,0.05,1,-45,0.05 


(Second  real  constant  set,  layers  3-4) 
R,  2 , 2 , ,  ,  ,  ,  , 

RMORE,  ,  ,  ,  ,  , 

RMORE, 2,-45,0.05,2,45,0.05 
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4.2  Results 

The  contour  plot  for  the  direction  of  the  axial  displacements  is  shown  in  Figure  3.  The  first 
item  to  note  is  how  the  displacements  behave  at  the  xz-symmetry  plane.  Instead  of  acting  like  the 
adjoining  row  of  elements,  the  displacements  increase  at  the  symmetry  plane.  This  behavior  should 
not  happen  at  a  symmetry  plane,  and  can  be  explained  by  the  fact  that  a  symmetry  plane  implies 
isotropy,  w'hich  is  not  the  case  for  lamina  with  fiber  orientations.  The  fiber  angles  for  each  layer  are 
mirrored  by  the  boundary  conditions  as  shown  in  Figure  4,  which  results  in  a  discontinuity.  At  this 
boundary  the  fiber  angles  are  reversed  which  is  why  the  displacements  are  reversed  in  reference  to 
the  free  edge  at  the  other  side  of  the  plate. 

The  second  feature  of  the  displacement  results  show  difference  in  displacements  between  the 
+45°  and  the  -45°  layers  through  the  plate  thickness.  This  can  be  seen  in  Figure  3.  and  Figure  5 
which  shows  the  contour  plot  for  displacements  perpendicular  to  the  load  direction.  This  is  expected 
behavior  for  the  given  laminate  ply  orientations. 

The  stress  output  was  not  incorporated  into  the  ANSYS  post-processor  due  to  time 
constraints.  The  stress  values  for  each  layer  is  available  in  tabuiar  form,  and  is  presented  in 
Appendix  B  for  the  third  column  of  elements  as  shown  in  Figure  2.  A  graph  of  the  stresses  across 
the  plate  width  for  the  second  lamina  layer  is  shown  in  Figure  6. 

The  results  of  Pipes  and  Pagano  revealed  that  at  the  free  edge  of  a  laminate  the  interlaminar 
stresses  (t„)  grow,  while  the  in-plane  shear  stresses  (rxy)  diminish  to  zero.  Also  there  is  a  decrease 
in  the  stresses  acting  along  the  load  direction  (o  j  at  the  free  edge.  The  results  shown  in  Figure  6 
follow  this  behavior,  but  with  some  discrepancies.  The  first  discrepancy  takes  place  from  the  xz- 
symrnetry  plane  to  y/b=0.3.  Here  the  stresses  are  effected  by  the  presence  of  the  symmetry  plane 
and  can  be  discounted.  The  second  discrepancy  occurs  near  the  free  edge,  in  that  -r,y  does  not 
diminish  to  zero,  end  that  xxz  does  not  increase  as  much  as  predicted  by  Pipes  and  Pagano.  This 
phenomena  can  be  explained  by  the  fact  that  an  8-noded  element  has  a  linear  displacement 
distribution  and  can  not  fully  capture  the  rapidly  changing  structural  response  at  the  free  edge.  Two 
methods  to  capture  this  response  is  to  use  a  finer  mesh  at  the  free  edge  and/or  use  20-noded  brick 
elements. 
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5.0  CONCLUSIONS 


The  Christensen  constitutive  equation  based  finite  element  method  has  successfully  been 
integrated  into  ANSYS  5.1.  The  results  follow  the  siress  distribution  at  the  laminate  free  edge 
according  to  the  results  of  Pipes  and  Pagano.  Therefore  with  proper  modeling  the  element  has  the 
accuracy  for  examining  potential  delaminations  at  design  features.  Proper  modeling  techniques 
include  not  using  symmetry  planes  and  using  a  fine  mesh  at  the  laminate  free  edge. 

The  problem  with  using  symmetry  planes  is  that  it  implies  isotropy,  which  was  shown  not 
to  be  the  case  for  the  element  formulation.  For  this  reason,  symmetry  planes  are  not  recommended 
for  modeling.  Further  development  of  the  element  formulation  algorithms  is  needed  to  take 
advantage  of  symmetry  planes. 

The  element  formulation  requires  that  there  are  two  material  sets  through  the  lamina 
thickness.  Using  only  one  material  set  gives  inaccurate  results.  Again  further  development  of  the 
element  formulation  algoritiims  are  needed  to  correct  this  deficiency. 

The  present  element  formulation  is  an  8  -noded  isoparametric  brick.  To  accurately  capture 
the  rapidly  changing  stress  field  near  the  laminate  free  edge,  the  mesh  will  have  to  be  finer  that  the 
rest  of  model.  Another  way  to  better  capture  this  stress  field  is  to  incorporate  a  20-noded  brick 
element  with  the  ANSYS  code. 

Writing  software  code  that  interfaces  with  an  existing  program  always  presents  challenges 
and  results  in  a  product  that  is  not  as  tightly  coupled  as  desired.  The  greatest  difficulty  in  writing 
the  user  element  was  dealing  with  the  input/output  with  ANSYS.  Determining  how  to  bring  in 
variables  from  ANSYS  and  stonng  results  was  not  intuitive  and  poorly  documented.  ANSYS  offers 
a  short  course  on  writing  user  programmable  features  and  in  retrospect  this  would  have  been  a 
valuable  course  for  developing  this  element  UPF. 
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6.0  RECOMMENDATIONS 


The  efficiency  of  using  this  element  could  be  improved  by  more  development  to  the 
algorithms  and  tighter  integration  with  ANSYS.  Specific  enhancements  are  allowing  the  user  to  use 
the  graphical  input  for  material  properties  and  real  constants,  adding  algorithms  to  take  advantage 
of  symmetry  planes,  refining  the  algorithms  so  that  only  one  material  property  is  needed, 
incorporating  20-node  brick  elements,  and  writing  the  code  to  the  view  the  strain  and  stress  output 
in  graphical  form. 

The  best  way  to  implement  these  enhancements  would  be  to  improve  the  element 
formulations  and  incorporate  them  into  the  existing  ANSYS  composite  element  SOLID  46.  This 
approach  would  give  the  tightest  coupling  between  the  Christensen  formulation  and  would  result  in 
the  most  efficient  use  of  this  element. 
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*DECK  UECI02  PARALLEL  USER  2/23/90  PCK 

SUBROUTINE  UECI02  (ELCDN.IELC.KERK) 

C  **"  THIS  SUBROUTINE  DEFINES  THE  CHARACTERISTICS  OF  US  I  RIO 
C  **♦♦  SEE  UEC100  FOR  DESCRIPTIVE  COMMENTS  ** 

C 

INCLUDE  TMPCOM' 

INCLUDE  ECHPRM 
C 

EXTERNAL  NMINFO,ALTINF,ANSERR,ERINQR,WRlNQR 
INTEGER  ER1NQR.WRINQR.IOTT 
INTEGER  IELC(*),I,KERR,KY2 
CHARACTER*  28  ELCDN 
C 

C  ***  ANSYS(R)  COPYRIGHT(O)  1971,78,82,83,85,87,89,92 

C  ***  SWANSON  ANALYSIS  SYSTEMS,  INC, 

C 

C  ****  DEFINE  ELEMENT  NAME  *** 

CALL  NMINFO  (IELC(1VUSER102 ') 

CALL  ALTINF  (IELC(  1 ),'  ') 

ELCDN  =  USER  ELEMENT  102 

****  ELEMENT  TATE  CHARACTERISTICS: 

****  BASED  ON  SOLID46  ELEMENT 

IELC(KDIM  )  =  3 
IELCdSHAP )  =  JBRICK 
IELC(IDEGEN)  =  3 
IELC(MSHLIM)  =  1 
IELC(KELSTO)  =  12 
IELC(KDOFS  )  =  UX  +  UY  +  UZ 
KY2  =  DELC(KYOP2 ) 

C 

IF  (KY2  .LE.  1)  THEN 
lELC(NMLAYS)  =  -1 
IF  (KY2  EQ.  1 )  ELC(NMLAYS)  =  -2 
IELC(LCANGL)  =  7 
lELC(NMDRLC)  =  312 
IELC(MATRQD)  =  1 
IELC(MATMUL)  =  1 
ELSE 

lELC(NMLAYS)  =  -3 
lELC(LCANGL)  =  128 
IELCCNMDRLC)  =  128 
END  IF 

IELC(NMNDMX)  =  8 

IELC(MATRXS)  =  STIFM  +  MASSM  +  SSTIFM 

IELC(NMPRES)  =  6 

IELC(NMIEMP)  =  8 

RETURN 

END 


SUBROUTINE  UEL102  (ELEM.IELC.ELMDAT.EOM  ASK. NODES, LOCSVR. KELREQ. 
XKELF1L.NR.XYZ.U.KELOUT.ZS.ZASS.DAMP.GST1F.ZSC.ZSCNR.ELVOL.ELMASS. 

X  CENTER, ELENER.ED1NDX.LCERST) 

C  —  SEE  EL  100  AND  EL101  FOR  DOCUMENTATION  — 

C 

C 

C  ***  ANSYS(R)  COPYRIGHT(C)  1971,78,82,83,85,87,89,92 

C  ***  SWANSON  ANALYSIS  SYSTEMS.  INC. 

C 

C  INPUT  ARGUMENTS; 

C  ELEM  (INT.SC.IN)  -  ELEMENT  LABEL  (NUMBER) 

C  IELC  (INT,AR(IELCSZ),IN)  -  ARRAY  OF  ELEMENT  TYPE  CHARACTERISTICS 
C  ELMDAT  (INT, AR(  1 0),IN)  -  ARRAY  OF  ELEMENT  DATA 

C  EOMASK  (INT, SC, IN)  -  BIT  PATTERN  FOR  ELEMENT  OUTPUT 

C  (SEE  OUTPCM) 

C  NODES  (INT, AR(NNOD), IN)  -  ARRAY  OF  ELEMENT  NODE  NUMBERS 
C  LOCSVR  (INT.SC.IN)  -  LOCATION  OF  THE  SAVED  VARIABLES 
C  ON  FILE  ESAV  FOR  THIS  ELEMENT 

C  KEI.REQCNl.AR(lO).IN)  -  MATRIX  AND  LOAD  VECTOR  FORM  REQUFSTS 
C  (INDICES  FOR  KELREQ  ARE  GIVEN  WITH  OUTPUT 

C  ARGUMENTS  BELOW) 

C  KELFIL  (INT,AR(10),IN)  -  KEYS  INDICATING  INCOMING  MATRICES  AND 

C  LOAD  VECTORS  (INDICES  FOR  KELFIL  ARE  THE 

C  SAME  AS  GIVEN  FOR  KELREQ  WITH  OUTPUT 

C  ARGUMENTS  BELOW) 

C  NR  (INT.SC.IN)  -  MATRIX  AND  LOAD  VECTOR  SIZE 

C  XYZ  (DP,AR(6,NNOD),IN)  -  NODAL  COORDINATES  (ORIG)  AND  ROTATION  ANGLE 
C  U  (DP,AR(NR,5),IN)  -  ELEMENT  NODAL  SOLUTION  VALUES 
C 

C  OUTPUT  ARGUMENTS: 

C  KELOUT  (INT,AR(  10), OUT)  -  KEYS  INDICATING  CREATED  MATRICES  AND 

C  LOAD  VECTORS  (INDICES  FOR  KELOUT 

C  ARE  THE  SAME  AS  FOR  KELREQ  BELOW) 

C  ZS  (DP. AR(NR, NRUNOUT)-  K  MATRIX  (KELREQ(l)) 

C  ZASS  (DP.AR(NR.NR).INOUT)-  M  MATRIX  (KELREQ(2)) 

C  DAMP  (DP,AR(NR,NR),INOUT)  C  MATRIX  (KJELREQ(3)) 

C  GSTIF  (DP,AR(NR,NR),lNOUT)-  S  MATRIX  (KELREQ(4)) 

C  ZSC  (DP,AR(NR),OUT)  -  APPLIED  F  VECTOR  (KELREQ(5)) 

C  ZSCNR  (DP.AR(NR).OUT)  •  N-R  RESTORING  F  VECTOR  (KELREQ(6)) 

C  OR  IMAGINARY  F  VECTOR  (KELREQ(7)) 

C  ELVOL  (DP, SC, OUT  -  ELEMENT  VOLUME 

C  ELMASS  (DP, SC, OUT)  -  ELEMENT  MASS 

C  CENTER  (DP, AR(3),OUT)  -  CENTROID  LOCATION 

C  ELENER  (DP,AR(5),OUT)  -  ELM  ENT  ENERGIES 

C  F.DINDX  (INT, AR(20),OUT)  -  ELEMENT  RES  ULT  DATA  FILE  INDEXES 

C  LCERST  (INT.SC.INOUT)  -  POSITION  ON  RESULT  FILE 

C 

C  -  START  OF  COMDECKS 
C 

INCLUDE  IMPCOM 
INCLUDE  ECHPRM 
INCLUDE  ELPARM 
INCLUDE  ELUCOM’ 

C 


C  — -  END  OF  COMDECKS 
C 

EXTERNAL  TRACK, ANSERR,SVGIDX,SVRGET.RVRGET,WRINQR, PROPEl, 
X VZERO 
C 

C  —  ARGUMENTS 

INTEGER  ELEM,IELC(IELCSZ),ELMDAT(10),EOMASK,NODES(8),LOCSVR, 
X  KELREQ(  10),  KELFIL{  1 0),NR,KELOUT(  10),EDINDX(2.0),LCERST.WRINQR 
C 

C  —  ARGUMENTS 
DOUBLE  PRECISION 

X  XYZ(6,8),U(NR,5)1ZS(NR,NR),ZASS(NK,NR),DAMP(NR1NR),GSTIF(NR,NR), 
XZSC(NR),ZSCNR(NR),ELVOL,ELMASS,CENTER(3)fLENER(5) 

C 

C  —  LOCAL 
C 

INTEGER 
X  SVTNDX(IO), 

XPMAT.IREAL.PREAL.NRVR.NSSVR.LLIOTT, 

X  K,KY2,NUMLAYERS,KPl,KP2 

C 

C 

C 

C  —  PROPERTIES 

C 

DOUBLE  PRECISION 
X  PROP(  10),PROPO, 

X  THK(IOO) 

C 

C  —  RVR.SSVR 
C 

DOUBLE  PRECISION 
X  RVR(312),  SSVR(10),SVIDNX(1),TK 
C 

INTEGER  ANGFIB(  1 00, 1 00),MATNUM( 250), THETA 
DOUBLE  PRECISION  PROPS(IO) 

C 

CALL  TRACK  (5.  UEL102') 

C 

C  —  DEFINE  INITIAL  DATA 
C 

MATNUM(ELEM)  =  ELMDAT(PMAT) 

IRF.AL  =  ELMDAT(PREAL) 

C 

C  —  IELC  POINTERS  DEFINED  IN  ECHPRM  AND  ELCCMT 
C 

NRVR  =  IELC(NMTRLC) 

NSSVR  =  IELC(NMSSVR) 

KY2  =  EELC(KYOP2) 

C 

C  DEFINE  KELOUT  SWITCHES  FOR  STIFFNESS  MATRIX 
C 

KELOUT(l)=l 

C 

C  DEFINE  KELREQ  SWITCHES  FOR  LOAD  VECTOR 


c 

KELREQ(5)  =1 
C 

C  —  GET  THE  S  VR  INDEX  VECTOR 
C 

CALL  SVGIDX  (LOCSVR.SVINDX) 

C 

C  —GET  THE  ELEMENT  REAL  CONSTANT  DATA 
C 

C  REAL  CONSTANT  DATA  IS  INPUTTED  USING  SAME  FORMAT 
C  AS  SOLID46  ELEMENT  WITH  KEYOPT(2)  =  0 
C  NL.LSYM.LP1  ,LP2,(BLANKS(2)),KREF,(BLANKS(5)), 

C  MAT, THETA, TK  FOB  LAYER  1 
C  MAT, THETA, TK  FOR  LAYER  2 
C  ETC.  UP  TO  LAYER  NL 
C 

C  FOR  THIS  USER  ELEMENT  HOWEVER,  THE  VARIABLES 
C  LSYM,LP1,LP2,BLANKS(2),KREF,BLANKS(5) 

C  WILL  NOT  BE  USED,  AND  WILL  INPUTTED  AS  BLANKS 

C 

C 

CALL  RVRGET  (ELEM,IREAL,IELC(  1  ),NRVR,RVR)C 
C 

IF  (KY2  ,EQ.  0)  THEN 
NUMLAYERS  =  RVR(l) 

K  =  13 
KP1  =  14 
KP2=  15 

DO  10 1  =  1,  NUMLAYERS 
MATNUM(I)  =  RVR(K) 

THETA  =  RVR(KPl) 

TK  =  RVR(KP2) 

C 

C  CONVERT  ANSYS  VARIABLES  TO  CHRISTENSEN  VARIABLES 
C 

ANGFIB(I,MATNUM(I))  =  THETA 
THK(I)  =  TK 
K  =  K  +  3 
KP1  =  KP1  +  3 
KP2  =  KP2  +  3 
10  CONTINUE 
ELSE 

WRITEdOTT.l) 

1  FORMAT  (/10X, '♦****  TAPERED  LAYER  OR  MATRIX  INPUT  NOT, 

X  /17X/YETLMPLEMENTED  ******  7) 

RETURN 
END  IF 
C 

C  —  SET  UP  MATERIAL  PROPERTIES 
r 

C  Ell 
C 

CALL  PROPE1  (ELEM.MATNUM(ELEM),  1 ,0.0,PROPO) 

PROPS() )  =  PROPO 
C 
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C  E22 
C 

CALL  PROPE 1  (ELEM  ,MATNUM(ELEM),2,0.0,PROPO) 

PR0PS(2)  =  PROPO 
C 

C  NU12 
C 

CALL  PROPE  1  (ELEM, MATNUM(ELEM), 4, 0.0, PROPO) 

PROPS(3)  =  PROPO 
C 

C  G12 
C 

CALL  PROPE  1(ELEM,MATNUM (ELEM), 7, 0.0, PROPO) 

PROPS(4)  =  PROPO 
C 

C  G13 
C 

CALL  PROPEl(ELEM,MATNUM(ELEM),9, 0.0,  PROPO) 

PROPS(5)  =  PROPO 
C 

C  DENSITY 
C 

CALL  PROPE  1  (ELEM, MATNUM (ELEM),  1 3,0.0, PROPO) 

PROPS(6)  =  PROPO 
C 

C  CHRIST.FOR  VARIABLE  P7  (LAMINA  THICKNESS)  IS  NOW 
C  INPUTTED  AS  A  REAL  CONSTANT  LN  ANSYS 
C  THEREFORE  PROPS(7)  WILL  BE  SKIPPED 
C  AND  SET  TO  ZERO 
C 

PROPS(7)  s  O.ODO 
C 

C  ALPHAX 

C 

CALL  PROPE1  (ELEM,  MATNUM(ELEM),  10, 0.0,  PROPO) 

PROPS(8)  =  PROPO 
C 

C  ALPHAY 
C 

CALL  PROPEl(ELEM,MATNUM(ELEM),l  1, 0.0, PROPO) 

PROPS(9)  =  PROPO 
C 

C  ALPHAZ 
C 

CALL  PROPEl(ELEM,MATNUM(ELEM), 12,0.0, PROPO) 

PROPS(10)  =  PROPO 
C 

C  CALCULATE  THE  STIFFNESS  MATRIX 
C 

CALL  BRICK  (ELEM,ZS,XYZ,NUMLAYERS,ANGFIB,PROPS,MATNUM) 

C 

C  WRITE  OUT  RESULTS 
C 

CALL  STRBRI(U,ELEM,XYZ,PROPS,NR,NODES,NUMLAYERS,MATNUM,ANGFIB, 
X  EDINDX.LCERST) 
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C 

CALL  TRACK  (15,  UEL102') 

RETURN 

END 

C 

Q  **************************************************************** 

c 

SUBROUTINE  BRICK(IELEM,EST1F, COORD, NLA YERS,IALPHA,PROPS,MATNO) 
C 

C  CALCULATES  THE  STIFFNESS  MATRIX  FROM  THE  INPUT  DATA 
C 

DOUBLE  PRECISION  ESTIF(24,24),BMATX(6,42),DMATX(6,6) 

DOUBLE  PRECISION  SHAPE(14),DERIV(3,14),DBMAT(6,42),CARTD(3,14) 
DOUBLE  PRECISION  ELCOD(6,8),COORD(6,8),TSTIF(33,33)  ,TEMP 
DOUBLE  PRECISION  POSGP(3),  WEIGP(3),  DVOLU 

COMMON  /CONTROL/  MPOIN,NPOIN,NELEM,NNODE,NDOFN,NDIMEtNSTRE, 

!  NPROP,NMATS,NVFIX,NEVAB,NINTR 

COMMON  /FIXEDBCV  PRESC(100,3),NOFIX(I00),IFPRE(100,3) 

DOUBLE  PRECISION  PROPS/ 10) 

COMMON  /MESH/  LNODS(6,8), 

!  ISNODE(750),ISELEM(500) 

EXTERNAL  TRACK,  WRINQR 

INTEGER  WRINQR.NLA  YERS  ,LALPHA(  1 00, 1 00),MATNO(250),1ELEM 
INTEGER  NGAUS 
CALL  TRACK  (6, 'BRICK') 

IOTT  =  WRINQR/2) 

C 

C 

C  EVALUATE  THE  COORDINATES  OF  THE  ELEMENT  NODAL  POINTS 
C 

C  —  SETS  UP  PARAMETERS  FOR  8-NODED  BRICK 
C 

NNODE  =  8 
NDIME  =  3 
C 

C  —  SETS  UP  CHRIST.FOR  ONLY  VARIABLES 
C 

NINTR=  1 
NEVAB  =  24 
NGAUS  =  3 
NSTRE  =  6 
NDOFN  =  3 
NPROP=  13 
C 

C  —  ZERO  OUT  VARIABLES 

C 

EXISP  =  0.0D0 
ETASP  =  0.0D0 
RHOSP  =  O.ODO 
DVOLU  =  O.ODO 
C 

C  —  SET  UP  GAUSSIAN  INTEGRATION  CONSTRANTS 
C 

CALL  GAUSSQ(NGAUS,POSGP,WEIGP) 

C 
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DO  10  INODE  =  l.NNODE 
DO  10  IDIME  =  l.NDIME 
ELCOD(IDIME, INODE)  =  COORD(IDIME.INODE) 

10  CONTINUE 
C 

C  INITIALIZE  THE  ELEMENT  STIFFNESS  MATRIX 
C 

DO  20  IEVAB  =  1  ,NEVAB+NINTR*9 
DO  20  JEVAB  =  1,NEVAB+NINTR*9 
TSTTF(  1EV  AB ,  JEV  AB )  =  0.0D0 
20  CONTINUE 
KGASP  =  0 
C 

C  ENTER  LOOPS  FOR  AREA  NUMERICAL  INTEGRATION,  STARTING  WITH  THE 
C  LOOP  FOR  EACH  LAYER  IN  THE  ELEMENT 
C 

DO  50 ILAYER  =  l.NLAYERS 
C 

C  EVALUATE  THE  D-MATRIX  IN  THIS  LAYER 
C 

CALL  DCHRIST(IELEM,ILAYER,IALPHA, PROPS, MATNO.DMATX) 

C 

C  NOW  THE  GAUSSIAN  INTEGRATION  OF  X  AND  Y  IN  THIS  LAYER 
C 

DO  50  IGAUS  =  1,NGAUS 
DO  50  JGAUS  =  l.NGAUS 
KGASP  =  KGASP+1 
EXISP  =  POSGP(IGAUS) 

ETASP  =  POSGP(iGAUS) 

RHOSP  =  - 1  .+(2.  ‘FLOAT (ILAYER)- 1  ,)/FLOAT(NLA  YERS) 

C 

C  EVALUATE  THE  SHAPE  FUNCTIONS,  ELEMENTAL  VOLUME,  ETC.  AT  THIS 
C  SAMPLING  POINT 
C 

CALL  SFBRIK(EXISP,ETASP,RHOSP,SHAPE,DERIV) 

CALL  JACBRIK(IELEM,DJACB,CARTD,DERIV,ELCOD) 

DVOLU  =  DJACB*WEIGP(IGAUS)*WEIGP(JGAUS) 

!  /FLOAT(NLAYERS) 

C  EVALUATE  THE  B  AND  BXD  MATRICES 
C 

CALL  BBRICK(BMATX.CARTD) 

CALL  DBE(BMATX,DMATX,DBMAT) 

C 

C  CALCULATE  THE  ELEMENT  STIFFNESSES 
C 

DO  30  IEVAB  =  1  ,NEVAB+NINTR*9 

DO  30  JEVAB  =  IEVAB, NEVAB+NINTR*9 
DO  30 ISTRE  =  l.NSTRE 
TSTIF(IEVAB, JEVAB)  =  TEMP  + 

!  BMATX(ISTRE,IEVAB)*DBMAT(ISTREdEVAB)*DVOLU 
30  CONTINUE 
50  CONTINUE 
C 

C  CONSTRUCT  THE  LOWER  TRIANGLE  OF  THE  STIFFNESS  MATRIX 
C 
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DO  60  IEVAB  =  1  ,NEVAB+NINTR*9 
DO  60  JEVAB  =  1,NEVAB+NINTR*9 
TSTIF(  JEVAB,  IEVAB)  =  TSTIF(IEVAB, JEVAB) 

60  CONTINUE 

C 

C  NOW  DO  THE  STATIC  REDUCTION  OF  TSTIF  TO  ESTIF 

C 

CALL  REDUCE  (TSTTF.ESHF.NINTR) 

C 

C  ALL  DONE 
C 

CAIX  TRACK(  16, ‘BRICK’) 

RETURN 

END 

C 

£  ************************************************************ 

C 

SUBROUTINE  REDUCE  (TSTIF, ESTIF, NINTR) 

C 

C  PERFORMS  A  STATIC  REDUCTION  OF  THE  TWO  INTERNAL  DEGREES  OF 

C  FREEDOM  IN  THE  8-NODE  BRICK  THAT  HAVE  BEEN  ADDED  TO  CAPTURE 
C  THE  THROUGH-THICKNESS  DISPLACEMENT  FIELD  IN  A  LAMINATED 

C  COMPOSITE  MATERIAL 
C 

EXTERNAL  WRINQR 
INTEGER  WRINQR 

DOUBLE  PRECISION  TSTIF(33,33),ESTIF(24,24),RSTIF(24,24),INDX(22) 
DOUBLE  PRECISION  STAU(22,24),STUA(24,22),STAA(22,22) 

DOUBLE  PRECISION  SINVAA(22,22),SINT(22,24) 

IOTT  =  WRINQR(2) 

C 

C  FIRST,  GET  THE  PART  OF  THE  MATRIX  THAT  HAS  JUST  THE 
C  ADDED  INTERNAL  DEGREES  OF  FREEDOM 
C 

NADD  =  9*NINTR 
C 

C  —  INITALIZE  ESTIF 
C 

DO  90  1=1,24 
DO  90  J=l,24 
90  ESTTF(I,J)  =  0.0D0 
C 

C  —  INITALIZE  STAA 
C 

DO  5  I  =  I, NADD 
DO  5  J  =  1,NADD 
STAA(I,J)  =  0.0D0 
5  CONTINUE 

IF  (NADD.GT.O)  THEN 

DO  10 IAA  =  I, NADD 

DO  10  KAA  =  l.NADD 

STAA(IAA.KAA)  =  TSTIF(IAA+24,KAA+24) 

10  CONTINUE 
C 

C  NOW  GET  THE  INVERSE  OF  THIS  MATRIX  FOR  THE  STATIC  REDUCTION 


,  • 


A- 8 


•  • 


I 

I 

« 


4 


« 

i 


i 

i 


i 


< 

l 


« 


C  EQUATION  -  NUMERICAL  RECIPES  IN  FORTRAN  -  P.  38 
C 

DO  20  IA  =  l.NADD 
DO  15  JA  =  l.NADD 

SINVAA(IA.JA)  =  0.0D0 
15  CONTINUE 

SINVAA(IA.IA)  =  1.0D0 
20  CONTINUE 

CALL  LUDCMP(STAA,NADD,22,INDX,D) 

DO  25  JA=  l.NADD 

CALL  LUBKSB(STAA,NADD,22,INDX,SINVAA(  1  ,JA)) 

25  CONTINUE 
C 

C  NOW  DO  THE  MATRIX  MUL  TIPLICATION  OF  [KUA]  X  [KAA  INVERSE]  X 
C 

C  FIRST  GET  [KUA]  AND  [KAU] 

C 

DO  30 IROW  =  l.NADD 
DO  30  ICOL  =  1,24 

STAU(IROW.ICOL)  =  TSHF(IROW+24,ICOL) 

30  CONTINUE 
DO  40  IROW  =  1.24 
DO  40  ICOL  =  l.NADD 

STUA(mOW.ICOL)  =  TSTIF(IROW,  ICOL-t-24) 

40  CONTINUE 
C 

C  NOW  MULTIPLY  OUT  THE  MATRICES 

C 

DO  50  mow  =  l.NADD 
DO  50 ICOL  =  1  24 
SINT(IROW.ICOL)  =  0.0D0 
DO  50  JROW  =  l.NADD 
SINTdROW.ICOL)  =  SINT(IROW,ICOL) 

!  +SINVAA(mOW,IROW)*STAU(mOW,ICOL) 

50  CONTINUE 

DO  60  mOW  =  1,24 
DO  60  ICOL  =  1,24 

RSHF(mow.icoL)  =  o.odo 

DO60  J  =  l.NADD 

RSTIFdROW.ICOL)  =  RSTDF(mOW.ICOL) 

!  +STUA(IROWJ)*$INT(J,ICOL) 

60  CONTINUE 
C 

C  NOW  SUBTRACT  RSTTF  FROM  TSTIF  TO  GET  ESTIF 

C 

DO  70  mow  =  1.24 
DO  70  ICOL  =  1,24 

ESTIF(IROW,ICOL)  =  TSTIF(mOW.ICOL)-RSTTF(mOW,ICOL) 

70  CONTINUE 
ELSE 

DO  80  mOW  =  1,24 
DO  80  ICOL  =  1.24 

ESTIF(mOW,ICOL)  =  TSTIF(mOW,ICOL) 

80  CONTINUE 
END  IF 


A-9 


c 

C  ALL  DONE 
C 

RETURN 

END 

C 

Q  ** a*********************************************  *v********* 

c 

SUBROUTINE  LUDCMP(A,N,NP,INDX,D) 

C 

C  ROUTINE  EXTRACTED  FROM  "NUMERICAL  RECIPES  IN  FORTRAN ",  PP.  35-36 
C 

PARAMETER  (NMAX=100,TIN  0=1.0-20) 

DOUBLE  PRECISION  A(NP,NP),INDX(NP),W(NMAX),AAMAX 
EXTERNAL  WRINQR 
INTEGER  WRINQR, NP 
IOTT  =  WRJNQR(2) 

D  =  l.ODO 
DO  12 1=1, N 

AAMAX  =  O.ODO 
DO  11  J=  1,N 

IF  (ABS(A(U)).GT.AAMAX)  AAMAX  =  ABS(A(U)) 

11  CONTINUE 

IF  (AAMAX.EQ.O.ODO)  THEN 

WRITE  (*,*) '  YOU  HAVE  A  SINGULAR  MATRIX" 

WRITE  (8,*)  A 
STOP 
ENDJF 

VV(I)  =  1. /AAMAX 

12  CONTINUE 
l,j  19  J  =  1,N 

DO  14  i  =  i,J-i 
SUM  =  A(I,J) 

DO  13  K  =1,1-1 

SUM  =  SUM  -  A(I,K)*A(K,J) 

13  CONTINUE 
A(I,J)  =  SUM 

14  CONTINUE 
AAMAX  =  O.ODO 

16  I=J,N 
SUM  =  A(U) 

DO  15  K  =  1,J-1 

SUM  =  SUM  -  A(I,K)*A(K,J) 

15  CONTINUE 
A(I,J)  =  SUM 

DrTI  =  VV(I)*ABS(SUM) 

B  DUM.GE.  AAMAX)  THEN 
IMAX  =  1 
AAMAX  -DUM 
ENDIF 

16  CONTINUE 

IF  (J.NE.EMAX)  THEN 
DO  17  K  =  1,N 

DUM  =  A(IMAX,K) 

A(IMAX.K)  =  A(J,K) 
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A(J.K)  =  DUM 

17  CONTINUE 
D  =  -D 

VV(IMAX)  =  VV(J) 

ENDIF 

INDX(J)  =  IMAX 

IF  (A(J,J).EQ,O.ODO)  A(J,J)  =  TINY 
IF  (J.NE.N)  THEN 
DUM  =  l./A(J.J) 

DO  18  I=J+1,N 

A(I,J)  =  A(I,J)*DUM 

18  CONTINUE 
ENDIF 

19  CONTINUE 
RETURN 
END 

C 

Q  ************************************************************ 

C 

SUBROUTINE  LUBKSB(A,N,NP,INDX,B) 

C 

C  EXTRACTED  FROM  "NUMERICAL  RECIPES  IN  FORTRAN"' 

C  DOES  THE  BACK  SUBSTITUTION  AFTER  THE  LU  DECOMPOSITION 
C 

DOUBLE  PRECISION  A(NP,NP),INDX(NP),B(NP) 

INTEGER  NP 
C 

11  =  0 

DO  121=  l.N 
LL  =  INDX(I) 

SUM  =  B(LL) 

B(LL)  =  B(I) 

IF  (II.NE.O)  THEN 
DO  11  J”  11,1-1 

SUM  =  SUM  -  A(I,J)*B(J) 

11  CONTINUE 

ELSEIF  (SUM.NE.0.0D0)  THEN 
11  =  I 

ENDIF 
B(I)  =  SUM 

12  CONTINUE 
DO  141  =  N.1,-1 

SUM  =  B(I) 

IF  (I.LT.N)  THEN 
DO  13  J  =  I+l.N 

SUM  =  SUM  -  A(I,J)*JB(J) 

13  CONTINUE 
ENDIF 

B(I)  =  SUM/A(I.I) 

14  CONTINUE 
RETURN 
END 

C 

Q  ************************************************************* 

c 
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SUBROUTINE  DBE(BMATX,DMATX,DBMAT) 

C 

C  CALCULATES  D  X  B 
C 

DOUBLE  PRECISION  BMATX(6.42),DMATX(6,6),DBMAT(6,42) 

COMMON  /control;  MPOIN.NPOIN.NELEM.NNODE.NDOFN.NDIME.NSTRE, 
!  NPROP,NMATS,NVFIX,NEVAB,NINTR 
C 

DO  10 ISTRE  =  l.NSTRE 
DO  10 IEVAB  =  1,NEVAB+9»NINTA 
DBMAT(ISTRE,IEVAB)  =  0.0D0 
DO  10  JSTRE  =  l.NSTRE 

DBMAT(ISTRE, IEVAB)  =  DBMAT(!STRE,IEVAB)+ 

!  DMATX(ISTRE,JSTRE)*EMATX(JSTRE, IEVAB) 

10  CONTINUE 
RETURN 
END 
C 

Q  I************************************************************* 

C 

SUBROUTINE  SFBRIK(S,T,R, SHAPE, DERTV) 

C 

C  CALCULATES  SHAPE  FUNCTIONS  AND  THEIR  DERIVATIVES  FOR 
C  BRICK  STRESS  AND  STRAIN  2-D  ELEMENTS 
C 

DOUBLE  PRECISION  SHAPE(14).DERIV(3,14) 

EXTERNAL  WRINQR 
INTEGER  WRINQR 
C 

C  —  INITIALIZE  SHAPE  FUNCTIONS  AND  DERIVATIVES 

C 

IOTT=WRINQR(2) 

DO  101=  1,14 
10  SHAPE(I)  =  0.0D0 
DO  201=  1.3 
DO  20  J=  1,14 
20  DERIV(I,J)=0.0D0 
C 

c 

SHAPE!  1)  =  0. 125*0  ,-T)*<  1  ,-S)*(l.-R) 

SHAPE(2)  =  0.125*0  ,-T)*(  1  .+S)*(  1  .-R) 

SHAPE(3)  =  0. 125*(  1  ,+T)*(  1  .+S)*(1.-R) 

SHAPE(4)  =  0. 1 25*(  1  ,+T)*(  1  .-S)*(  1  .-R) 

SHAPE(5)  =  0. 1 25  *(  1  ,-T)*(  1  ,-S)*(  1  ,+R) 

SHArL(6)  =  0. 1 25*(  1  ,-T)*(  1  ,+S)*(  1  .+R) 

SHAPE(7)  =  0. 1 25*(  1  ,+T)*(  1  ,+3)*(  1  ,+R) 

SHAPEC8)  =  0. 1 25*(  1  ,+T)*(  1  ,-S)*Q  .+R) 

SHAPE(9)  =  l.-S**2 
SH.\TE(10)=  l.-T**2 
SHAPE(ll)  =  l.-R**2 
SHAPE(12)  =  S-S**3 
SHAPE(13)  =  T-T**3 
SHAPE(14)  =  R-R**3 
C 


A-12 


on  non  non 


( 


AND  THEIR  DERIVATIVES 

DERIV(l.I)  =  -,125*(1.-T)*(1.-R) 
DERIV(1,2)  =  125*(1.-T)*(1  -R) 
DERIV(l.j)  =  .125*(1.+T)M.-R) 
DERIV(1,4)  =  -.125*0.+T)M.-R) 
DERJV(1,5)  =  -.125*(1.-T)*(1.+R/ 
DERTV(1,6)  =  ,125*(1.-T)*(1.+R) 
DERIV;i,7)  =  ,125‘'(1.+T)*(1.+R) 
DERTV(1,8)  =  -.125*(1.+T)*(1.+R) 
DER1V(1,9)  =  -2.*S 
DERIV(1,10)  =  O.ODO 
DERIV(1,11)  =  0.0D0 
DERIV(1,12)  =  l.-3.*S**2 
DERIV(1,13)  =  0.0D0 
DERIV(1,14)  =  0.0D0 
DERIV(2,1)  =  -125*(1.-S)*(1.-R) 
DERIV(2,2)  =  -,125*(1.+S)*(1.-R) 
DERIV(2,3)  =  ,125*(1.+S)*(1.-R) 
DERIV(2.4)  =  ,125*(1.-S'*(1.-R) 
DERTV(2,5)  =  -.125*(1.-S)*(l.+R) 
DERIV(2.6)  =  -.  1 25  ♦(  1  ,+S)*(  1  ,+R) 
DERIV(2,7)  =  ,125*(1.+S)*(1.+R) 
DERIV(2.8)  =  .125*(1.-S)*(1.+R) 
DERTV<2,9)  =  0.0D0 
DERJV(2,10)  =  -2.*T 
DERTV(2,1 1)  =  O.ODO 
DERTV(2,12)  =  O.ODO 
DERTV(2,13)  =  l.-3.*T**2 
DERW2,14)  =  O.ODO 
DERIV(3,1)  =  -.125*(1.-S)*(1.-T) 
DERIV(3,2)  =  -.125*(1.+S)*(!.-T) 
DERTV(3,3)  =  -,125*(1.+S)*(1.+T) 
DERIV(3,4)  =  -.125*(1.-S)*(1.+T) 
DERIV(3.5)=  .125*(1.-S)*(1.-T) 
DERIV(3,6)  =  .125*(1.+S)*(1.-T) 
DERTV(3,7)  =  ,125*(1.+S)*(1.+T) 
DERIV(3,8)  =  .125*(1.-S)*(1.+T) 
DERTV(3,9)  =  O.ODO 
DERIV(3.I0)  =  O.ODO 
DERIV(3,1 1)  =  -2.*R 
DERIV(3,12)  =  O.ODO 
DERIV(3,13)  =  O.ODO 
DERIV(3,14)  =  l.-3.*R**2 

ALL  DONE 

RETURN 

END 
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SUBROUTINE  JACBRIK(IELEM,DJACB.CARTD,DERIV(ELCOD) 

CALCULATES  THE  JACOBIAN  MATRIX  AND  DETERMINANT  AND  INVERSE 
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C  FOR  3-D  ELEMENTS 
C 

DOUBLE  PRECISION  XM(3.3).XI(3,3),CARTDi3,14) 

C 

COMMON  /CONTROL/  MPOIN.NPOIN,NELEM,NNODE,NDOFN,NDIME,NSTRE, 
!  NPROP.NMATS.NVFIX.NEVAB.NINTR 

DOUBLE  PRECISION  DERIV(3,14).ELCOD(6,8) 

INTEGER  IELEM 
EXTERNAL  WRINQR 
INTEGER  WRINQR 
IOTT  =  WRINQR(2) 

C 

C  —  INITIALIZE  XI  MATRIX 
C 

DO  101=  l.NDIME 
DO  10  J=  l.NDIME 
XM(I,J)  =  0.0D0 
10  XI(I.J)  =  0.0D0 
D.IACB  =  0.0D0 
C 

C  CREATE  JACOBIAN  MATRIX 
C 

DO  20 IDIME  =  l.NDIME 
DO  20  JDIME  =  l.NDIME 
DO  20  INODE  =  1  .NNODE 

XM(IDIME.JDIME)  =  XM(IDIME.JDIMEj+ 

!  DERTV(IDIME,INODE)*ELCOD(JDiME.INODE) 

20  CONTINUE 
C 

C  NOW  THE-  DETERMINANT  AND  INVERSE  OF  THE  JACOBEAN 
C 

DJACB  =  XM(1,1)*(XM(2,2)*XM(3,3)-XM(3,2)*XM(2,3)) 

!  -XM(  1 ,2)*(XM(2, 1  )*XM(3.3)-XM(3, )*XM(2,3)) 

!  +XM(  1 ,3)*<XM(2, 1  )*XM(3,2)-XM(3, 1  )*XM(2,2)) 

IF  (DJACB.GT.0.)  GO  TO  30 
WRITE  (IOTT, 999)  IELEM 

999  FORMAT  (//.’  ELEMENT  NUMBER  .15.  HAS  ZERO  OR 
!  NEGATIVE  AREA  -  PLEASE  CHECK  INPUT////) 

STOP 

30  XI(l.l)  =  (XM(2,2)*XM(3,3)-XM(3,2)*XM^2,3))/DJACB 
XI(  1,2)  =  (XM(3,2)*XM(1,3)-XM(1,2)*XM(3,3))/DJACB 
XI(  1 ,3)  =  (XM(  1 ,2)*XM(2,3)-XM(2.2)*XM(  1 ,3))/D  j  ACB 
XJ(?.,1)  =  (XM(3, 1)*XM(2,3)-XM(2, 1  )*XM(3,3))/DJACB 
XI(2,2)  =  (XM(  1 , 1)*XM(3,3)-XM{3, 1  )*XM(  1 ,3))/DJACB 
XI(2,3)  =  (XM(2,l)*XM(1.3)-XM(l,l)*XM(2,3))/DJACB 
XI(3,1) « (XM(2.1)*XM(3,2)-XM(3,1)*XM(2,2)VDJACB 
XI(3.2)  *  fXM(3,l)*XM(l,2)-XM(l,l)*XM(3,2))/DJACB 
XI(3,3)  =  (XM(  1 , 1  )*XM(2,2)-XM(2, 1  )*XM(  1 ,2)VDJ  ACB 
C 

C  CALCULATE  CARTESIAN  DERIVATIVES 
C 


50  CARTD(I.J)  =  O.ODO 
DO  40  ID1ME  =  l.NDIME 
DO  40  INODE  =  l.NNODE+NINTR*3 
C  CARTD(IDIME.INODE)  =  O.ODO 
DO  40  JDIME  =  l.NDIME 

CARTDdDIME.INODE;  =  CARTD(IDIME.1N0DE)+ 

!  XI(.  JlME.JDIME)*DERIV(JDIME,INODE) 

40  CONTINUE 
C 

C  ALL  DONE 
C 

RETURN 

END 

C 

q  ****************** ******** ************* ***************** 

C 

SUBROUTINE  BBRICK(BMATX.CARTD) 

C 

C  CALCULATES  THE  STRAIN  MATRIX  B  FOR  THE  BRICK  ELEMENT 
C 

DOUBLE  PRECISION  BMATX(6.42),CARTD(3.14) 

COMMON  /CONTROL/  MPOIN.NPOIN.NELEM.NNODE.NDOFN.NDIME.NSTRE, 
!  NPROP.NMATS.NVFIX.NEVAB.NINTR 
C 

DO  5  1=1,6 
DO  5  J  =  1,42 
5  BMATX(I.J)  =  O.ODO 

DO  10  INODE  =  l,NNODE+NINTR*3 

MGASH  =  (INODE-l)*3+l 

NGASH  =  (INODE- 1)*3+2 

LGASH  =  (INODE- 1)*3+3 

BMATX( ! , MGASH)  =  CARTD'  1, INODE) 

BMATX(l.NGASH)  =  O.ODO 
BMATX(l.LGASH)  =  O.ODO 
BMATX(2, MGASH)  =  O.ODO 
BMATX(2,NGASH)  =  CARTD(2  INODE) 

BMATX(2,LGASH)  =  O-ODO 
BMATX(3, MGASH)  =  O.ODO 
BMATX(3,NGASH)  =  O.ODO 
BMATX(3,LGASH)  =  CARTD(3,INODE) 

BMATX(4, MGASH)  =  CARTD(2.INODE) 

BMATX(4,NGASH)  =  CARTD(1,  INODE) 

BMATX(4,LGASH)  =  O.ODO 
BMATX(5, MGASH)  =  O.ODO 
BMATX(S.NGASH)  =  C  ARTD(  3  .INODE) 

BMATX(5,LGASH)  =  CARTD(2,INODE) 

BMATX(6,MGASH)  =  CARTD(3,INODE) 

BMATX(6,NGASH)  =  O.ODO 
BMATX(6,LGASH)  =  CARTD(1, INODE) 

10  CONTINUE 
C 

C  ALL  DONE 
C 

RETURN 

END 
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c 

SUBROUTINE  DCHRIST(IELEM.ILAYER.IALPHA.PROPS.MATNO.DMATX) 

c 

C  ROUTINE  TO  EVALUATE  ELASTICITY  MATRIX  FOR  THE  ITH  LAYER  OF 
C  THE  ELEMENT.  ACCORDING  TO  CHRISTENSEN  ( 1 990)  EQUATION  25 

C 

COMMON  /CONTROL/  MPOIN.NPOIN.NELEM.NNODE.NDOFN.NDIME.NSTRE, 

!  NPROP.NMATS.NVFIX.NEVAB.NINTR 

COMMON  /FIXEDBC7  PRESC(100,3),NOFIX(100),IFPRE(100,3) 

INTEGER  LALPHA(  100, 100),MATNO( 250),IELEM 

DOUBLE  PRECISION  PROPS(10).AVEMU.AVELAM,AVEMOD,A1,A2 

COMMON  /MESH/  LNODS(6,8), 

!  ISNODE(750),ISELEM(500) 

C 

DOUBLE  PRECISION  DMATX(6,6) 

C 

C 

C  MATERIAL  PROPERTIES  INPUT  AS  FOLLOWS 
C 

C  PI  =  El  1 
C  P2  =  E22 
C  P3  =  NU 1 2 
C  P4  =  G12 
C  P5  =  G13 
C  P6  =  DENSITY 
C  P7  =  LAMINA  THICKNESS 
C  P8  =  THERMAL  EXPANSION  -  X 
C  P9  =  THERMAL  EXPANSION  -  Y 
C  P10  =  THERMAL  EXPANSION  -  Z 

C 

C  CALCULATE  THE  MATRIX  DOMINATED  PROPERTIES  AND  THE  AXIAL  PROPERTY 

C 

C 

C  —  ZERO  OUT  VARIABLES 
C 

AVEMU  =  O.ODO 
AVELAM  =  O.ODO 
AVEMOD  =  O.ODO 
A1  =  O.ODO 
A2  =  O.ODO 
C 

IF  (PROPS(5).LE.l.)  THEN 

AVEMU  s  (PROPS(2)*(1.-PROPS(3))/(2.*(1.-PROPS(2)/ 

!  PROP$(l)*PROPS(3)**2))+PROPS(4))/2, 

ELSE 

AVEMU  =  (PROPS(2)*(  1  .-PROPS(3))/(  1  ,-PROPS(2)/PROPS(  1  )* 

!  PROPS(3)**2)+2.+PROPS(4)+PROPS(5))/5. 

ENDIF 

AVELAM  =  AVEMU*2.*PROPS(3)/(l.-2.*PROPS(3)) 

AVEMOD  =  PROPS(  l)-AVEMU*2.*(l.+PROPS(3)) 

C 

C  ZERO  OUT  THE  D-MATPiX 
C 
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DO  10  ISTRE  =  l.NSTRE 
DO  10  JSTRE=  l.NSTRE 
DMATX(ISTRE.JSTRE)  =  0.0D0 
10  CONTINUE 

EVALUATE  THE  DIRECTION  COSINES  FOR  THIS  LAYER 
CONVER  =  3.1415926536/180. 

A1  =  COS(CONVER*FLOAT(IALPHA(ILAYER.MATNO(IELEM)))) 

A2  =  SIN(CONVER*FLOAT(IALPHA(ILAYER,MATNO(IELEM)))) 

NOW  THE  COMPONENTS  OF  THE  D-MATRIX  FOR  THIS  LAYER 

DMaTX(I.I)  =  AVELAM*(l.-PROPS(3))/PROPS(3)  +  AVEMOD*Al**4 
DMATXG.2)  =  AVELAM  +  AVEM0D*A1**2*A2»*2 
DMATX0.3)  =  AVELAM 
DMATX0.4)  =  AVEM0D*A1**3*A2 
DMATX(2,1)  =  DMATX(1,2) 

DMATX(2,2)  =  AVELAM*(1.-PROPS(3))/PROPS^)  +  AVEMOD*A2**4 
DMATX(2.3)  =  AVELAM 
DMATX(2,4)  =  AVEMOD*Al*A2**3 
DMATX(3,1)  =  DM  ATX(  1 ,3) 

DMATX(3,2)  =  DMATX(2,3) 

DMATX(3.3)  =  AVELAM*(l.-PROPS(3))/PROPS(3) 

DMATX(4,1)  =  DMATX(1.4) 

DMATX(4,2) «  DMATX(2,4) 

DMATX(4,4)  =  AVEMU  +  AVEMOD*Al**2*A2**2 
DMATX(5,5)  =  AVEMU 
DMATX(6,6)  =  AVEMU 

ALL  DONE 

RETURN 
END 

**•  *********************  I************************************ 

SUBROUTINE  GAUSSQ(NGAUS,POSGP,WEIGP) 

ROUTINE  TO  SET  UP  SAMPLING  POINTS  AND  WEIGHTING  FACTORS 
FOR  THE  ELEMENT  NUMERICAL  INTEGRATIONS 

COMMON  /CONTROL/  MPOIN,NPOIN,NELEM,NNODE,NDOFN,NDIME,NSTKE, 
!  NPROP.NMATS.NVFIX,NEVAB,NINTR 
DOUBLE  PRECISION  POSGP(?),WEIGP(3) 

INTEGER  KGAUS.IGAUS.JGAUS 
C 

EXTERNAL  TRACK,  WRINQR 
INTEGER  WRINQR,  NGAUS 
IOTT  =  WRINQR(2) 

C 

CALL  TRACK(7,’GAUSSQ) 

IF  (NGAUS.GT.2)  GO  TO  10 
POSGP(l)  =  -.577350269189626 
WEIGP(')  =  1.0DC 
GO  TO  20 


A-17 


nnnnn 


10  POSGPU)  =  . 774596669241483 
C 

C  CHANGE  BY  M.  JONES  ON  1/15/96 

C  CHANGED  POSGP(3)  TO  POSGP(2)  AND  WEIGP(3)  TO  WEIGP(2) 
C  TO  REFLECT  GUASS-QUAD  SAMPLING  POINTS  AND  WEIGHTS 
C  WHEN  KGUAS  =  3 
C 

PO$GP(2)  =  O.ODO 
WEIGP(l)  =  .555555555555556 
WEIGP(2)  =  O.ODO 
20  KGAUS  =  NGAUS/2 
DO  30 IGASH  =  l.KGAUS 
JGASH  =  NGAU  S+ 1  -IGASH 
POSGP(JGASH)  =  -POSGP(IGASH) 

WEIGP(JGASH)  =  WEIGP(IGASH) 

30  CONTINUE 
CALL  TRACK!  17,'GAUSSQ) 

RETURN 

END 


Mi********************************************************* 

SUBROUT J'JESTRBRI(ASDIS,IELEM, COORD, PROPS,NR,LNODS,NLAYERS, 

X  MATNO.IALPHA.LCERST.EDINDX) 

C 

C  ROUTINE  TO  CALCULATE  STRESSES  IN  BRICK  ELEMENTS  AFTER 
C  DISPLACEMENTS  HAVE  BEEN  CALCULATED 
C 

C  LOGICAL  THERMAL 

C  DOUBLE  PRECISION  STRIN(6,8),STRAN(6,8) 

EXTERNAL  WRINQR,UEP102,TRACK,ELDWRT 

INTEGER  WRINQR.NLA  YERS  ,IALPHA(  1 00, 1 00),MATNO(250).IELEM,NR 
INTEGER  NGAUS.SVINDX(l) 

DOUBLE  PRECISION  ELDIS(3,8),STRES(6),STEMP(27,6),STRAI(6), 

X  SRTMP(27,6),SMATX(6,24,27),SRMAT(6,24,27),BMATX(6,42), 

X  DMATX(6,6) 

DOUBLE  PRECISION  DBMAT(6,42),CARTD(3,14),COORD(6,8),ELCOD(6,8; 
DOUBLE  PRECISION  SHAPE(14),DERIV(3,14),PROPS(10) 

DOUBLE  PRECISION  ASDIS(NR,5) 

DOUBLE  PRECISION  POSGP(3),WEIGP<3) 

INTEGER  EDINDX(20),LCERST 
INTEGER  LNODS(8) 

INCLUDE  ELPARM' 

C 

COMMON  /CONTROL/  MPOIN,NPOIN,NELEM,NNODE,NDOFN,NDIME,NSTRE, 

!  NPROP,NMATS,NVFEX,NEVAB,NINTR 

COMMON  /FIXEDBC/  PRESC(100,3),NOFDC(100),IFPRE(100,3), 

!  ISNODE(750),ISELEM(500) 

C  KGAST  =  0 
C  THERMAL  =  FALSE. 

C 

NNODE  =  8 
NDIME  =  3 
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NINTR  =  l 
NEVAB  =  24 
NGAUS  =  3 
NSTRE  =  6 
NDOFN  =  3 
NPROP  =  13 
C  WRITE  HEADER 
C 

IOTT  =  WRINQR(2) 

WRITE  (IOTT,*)  •  CALCULATING  STRESSES' 

WRITE  (6,2) 

2  FORMAT  (/,'  ELEM.2X, 

!  '  X  Y  Z'j:  SX  SY  SZ'. 

!  1  SXY  SYZ  SXZ) 

c 

C  CALCULATE  AND  OUTPUT  STRESSES  FOR  EACH  ELEMENT 
C 

KGAST  =  0 
C 

C  EVALUATE  THE  COORDINATES  OF  THE  ELEMENT  NODES 
C 
C 
C 

DO  80  INODE  =  l.NNODE 
C  LNODE  =  IABS(LNODS(IELEM,  INODE)) 

DO  80 IDIME  =  l.NDIME 

ELCODdDIME, INODE)  =  COORD(IDIME, INODE) 

80  CONTINUE 
XCENT  *  O.ODO 
YCENT  =  O.ODO 
C 

C  —  Z-HEIGHT  OF  ELEMENT 
C 

DZELEM  >  ELCOD(3,5  )-ELCOD(3, 1 ) 

Z1  =ELCOD(3.I) 

DO  9  INODE  =  l.NNODE 

XCENT  =  XCENT  +  ELCOD(l,INODE)M25 
YCENT  =  YCENT  +  ELCOD(2,INODE)M25 

9  CONTINUE 
C 

C  IDENTIFY  THE  DISPI ACEMENTS  OF  THE  ELEMENT  NODAL  POINTS 
C 

NPOSN  =  0 

DO  10  INODE  =  l.NNODE 
DO  10 IDOFN  =  1 .NDOFN 
NPOSN  =  NPOSN+1 

ELDIS(IDOFN, INODE)  =  ASDIS(NPOSN,l) 

10  CONTINUE 
C 

C  CALCULATE  THE  STRESS  AND  STRAIN  MATRIX  FOR  THE  ELEMENT 
C 

DO  90  ILAYER  =  l.NLAYERS 
RHOSP  =  -l.+{2  *FLOAT(ILAYER)-1.)/FLOAT(NLAYERS) 

ZCFNT  =  Z1  +  .5*(l.+RHOSP)*DZELEM 
C 
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C  ZERO  OUT  THE  TEMPORARY  STRESSES 
C 

NGASP  =  NGAUS’NGAUS 
DO  5  JSTRE  =  l.NSTRE 
DO  5  IGASP  =  1 , NGASP 
STEMP(IGASP.JSTRE)  =  0. 

SRTMP(IGASP.JSTRE)  =  0. 

5  CONTINUE 
C 

C  EVALUATE  THE  D-MATRIX 
C 

CALL  DCHRIST(IELEM, ILAYER.IALPHA, PROPS  .MATNO.DMATX) 

C 

KGASP  -  0 

DO  40 IGAUS  =  1  .NGAUS 
DO  40  JGAUS  =  1  NGAUS 
KGASP  =  KGASP+I 
EX3SP  =  POSGP(IGAUS) 

ETASP  *  POSGP(JGAUS) 

C 

C  EVALUATE  THE  SHAPE  FUNCTIONS,  ELEMENTAL  VOLUME,  ETC. 

C 

CALL  SFBRIK(EXISP,ETASP,RHOSP,SHAPE,DERIV) 

CALL  JACBRIK(IELEM.DJACB,CARTD.DERIV1ELCOD) 

C 

C  EVALUATE  THE  B  AND  BXD  MATRICES 
C 

CALL  BBRICK(BMATX.CARTD) 

CALL  DBE(BMATX,DMATX,DBMAT) 

C 

C  STORE  THE  COMPONENTS  OF  THE  STRESS  AND  STRAIN  MARTICES 
C 

DO  15  ISTRE  =  l.NSTRE 
DO  15  IEVAB  =  l.NEVAB 

SMATX(ISTRE,IEVAB, KGASP)  =  DBMAT(ISTRE, IEVAB) 

SRMAT(ISTRE, IEVAB, KGASP)  =  BM/,TX(IS1REJEVAB) 

15  CONTINUE 
40  CONTINUE 
C 

C  NOW  TO  CALCULATE  STRESSES 
C 

KGASP  =  0 
C 

C  ENTER  LOOPS  OVER  EACH  GAUSS  POINT  IN  THE  LAYER 
C 

DO  50  IGAUS  =  1, NGAUS 
DO  50  JGAUS  =  1  NGAUS 
KG  AST  =  KGAST+1 
KGASP  =  KGASP+1 
C 

C  COMPUTE  THE  STRESS  AND  STRAIN  COMPONENTS  AT  THE  SAMPLING 
C  POINTS 
C 

DO  20 ISTRE  =  l.NSTRE 
KGASH  =  0 
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DO  20  INODE  =  l.NNODE 
DO  20  IDOFN  =  1  .NDOFN 
KG  ASH  =  KGASH+1 

STEMP(KGASP.ISTRE)  =  STEMP(KGASP,ISTRE)+ 

!  SMATX(ISTRE,KGASH,KGASP)*ELDIS(IDOFN, INODE) 
SRTMP(KGASP.ISTRE)  =  SRTMP(KGASP,ISTkE)+ 

!  SRMAT(ISTRE,KGASH,KGASP)*ELDIS(IDOFN, INODE) 

20  CONTINUE 
C 

C  FOR  THERMAL  LOADING,  ADD  ON  THE  INITIAL  THERMAL  STRESS 
C 

C  IF  (THERMAL)  THEN 
C  READ  ( 1  ,REC=IELEM)  STRIN.STRAN 
C  DO  30  IS'TRI  =  l.NSTRE 

C  STEMP(KGASP.ISTRl)  =  STEMP(KGASP,ISTR1)+STRIN(ISTR1,KGAST) 
C  SRTMP(KGASP.ISTRl)  =  SRTMP(KGASP,ISTRI)+STRAN(ISTR1,KGAST) 
C  30  CONTINUE 
C  END  IF 
50  CONTINUE 

3  FORMAT  (214, 3F10.4,2(/,1X,6G12.5)) 

OUTPUT  STRESSES  FOR  THIS  LAYER 


—  AT  THIS  POINT  UEP102  SHOULD  BE  CALLED  FOR  OUTPUT 

NGASP  =  NGAUS*NGAUS 

DGASP  =  2.  *FLO  AT(NGASP) 

DO  60  LSTRE=I,NSTRE 
STRES(LSTKE)  =  0. 

STRAI(LSTRE)  =  0. 

60  CONTINUE 
DO  70  KSTRE  =  l.NSTRE 
DO  70 IGASP  =  1, NGASP 

STRES(KSTRE)  =  STRES(KSTRE)+STEMP(IGASP,KSTRE)/DGASP 
STRAI(KSTRE)  =  STRAI(KSTRE)+SRTMP  (IGASP, KSTRE)/DGASP 
70  CONTINUE 

90  WRITE  (IOTT.3)  IELEM,ILAYER,XCENT,YCENT,ZCENT,(STRES(ISTRE), 
!  ISTRE=  1  ,NSTRE),(STRAI(ISTRE),ISTRE=  1  ,NSTRE) 

CALL  ELDWRTdELM, EDENS, LCERST,EDINDX(1),NSTRE, SIRES) 

RETURN 

END 


First  Row  of  Elements  (45,-45  layers) 


ELMDAT(  1)  =  1  for  Element  45 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 

45  1  0.5000  -0.0500  0.0250 
85665.  -487.81  481.11  24702.  82.856  8328.7 

0.3095 1E-0 1 -0.20 1 27E-0 1 -0.26677E-02-0.44800E-02  0.98246E-04  0.98758E-02 
45  2  0.5000  -0.0500  0.0750 

87187.  1504.4  521.47  -26619.  -69.370  8306.2 

0.30844E-01-0.19955E-01-0.26677E-02  0.41761E-02-0.82256E-04  0.98491E-02 

ELMDAT(  1)=  1  for  Element  39 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 

39  1  0.5000  -0.1500  0.0250 
91559.  687.80  -355.99  35803.  49.614  563.02 

0.30953E-01-0.22922E-01-0.22899E-02-0.48002E-04  0.58830E-04  0.66760E-03 
39  2  0.5000  -0.1500  0.0750 
Ronn  .0003 «  4m  .iisn 

0.30985E-01-0.23030E-01-0.22899E-02  0.54083E-03-0.13846E-03  0.66295E-03 


ELMDAT(  1)  =  1  for  Element  33 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 
33  1  0.5000  -0.2500  0.0250 
91228.  -503.58  43.513  34929.  -43.803  50.369 

0.31 176E-01-0.23210E-01-G.20985E-02-0. 15644E-03-0.51939E-04  0.59726E-04 
33  2  0.5000  -0.2500  0.0750 
92541.  511.04  46.429  -36300.  -38.825  45.634 

0.31266E-01-0.23296E-01-0.209S5E-02-0.95467E-04-0.46037E-04  0.541 1  IE  04 


ELMDAT(  1)=  1  for  Element  27 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 

27  1  0.5000  -0.3500  0.0250 

92021.  192.20  60.987  35598.  15.948  -0.80817 

0.3 1 249E-0 1  -0.23 195E-0 1  -0.2 1 144E-02-0. 1054 1 E-03  0. 1 89 1 1 E-04-0.95829E-06 

27  2  0.5000  -0.3500  0.0750 

93056.  1234.2  69.668  -36795.  24.362  -10.832 

0.3 1 254E-0 1  -0.23 1 86E-0 1  -0.2 1 144E-02-0. 1 0709E-03  0.28887E-04-0. 1 2844E-04 

ELMDAT(  1)=  1  for  Element  21 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 

21  1  0.5000  -0.4500  0.0250 

92654.  999.15  69.529  36388.  18.344  -1.2282 

0.3 1 210E-0 1 -0.23 1 30E-0 1-0.2 1 176E-02  0.20634E-04  0.2 1 75  IE-04-0. 14564E-05 

21  2  0.5000  -0.4500  0.0750 

92285.  732.56  48.222  -36100  12.680  -12.178 

0.31 162E-01-0.231 17E-01-0.21 176E-02  0.401 14E-05  0.15035E-04-0.14440E-04 

ELMDAT(  1)=  1  for  Element  15 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 

15  1  0.5000  -0.5500  0.0250 

92654.  1148.9  82.563  36468.  27.817  -2.3970 

0.31 1 67E -0 1 -0.23085E-01  -0.2 1 125E-02  0.33844E-04  0.32984E-04-0.28422E-05 

15  2  0.5000  -0.5500  0.0750 

91328.  -39.729  -12.260  -35374.  29.173  -9.7503 

0.31 048E-0 1  -0.23 1 2 IE-01-0.2 1 1 25E-02  0.40684E-04  0. 34592E-04-0.il 56 IE-04 


ELMDAT(  1)  =  1  for  Element  9 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 

9  1  0.5000  -0.6500  0.0250 

91468.  411.51  -311.82  35112.  -172.14  -443.69 

0.31 105E-01-0.22880E-01-0.23220E-02-0.34010E-03-0.20412E-03-0.52610E-03 
9  2  0.5000  -0.6500  0.0750 

90467.  -266.98  -471.66  -34788.  3.1402  -438.84 

0.30879E-0 1 -0.229 1 6E-0 1  -0. 23220E-C2  0.18055E-03  0.37235E-05-0.52036E-03 

ELMDAT(  1)=  1  for  Element  3 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 

3  1  0.5000  -0.7500  0.0250 

86481.  279.66  250.74  25474.  -189.82  -8152.0 

0.3101 6E-0 1 -0,2009 1 E-0 1 -0.27950E-02-0.44209E-02-0. 22508E-03-0.96663E-02 
3  2  0.5000  -0.7500  0.0750 

85291.  -459.24  51.317  -25167.  33.788  -8122.8 

0.307 1 9E-0 1 -0.20 1 20E-0 1 -0.27950E-02  0.42034E-02  0.40064E-04-0.96317E-02 


Second  Row  of  Elements  (-45,45  layers) 


ELMDAT(  2)  =  2  for  Element  93 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 
93  1  0.5000  -0.0500  0.1250 

87187.  1504.4  521.47  -26619.  69.370  -8306.2 

0.30844E-0 1  -0. 1 9955E-0 1-0.26677E-02  0.4 1 76 1 E-02  0. 82256E-04-0.9849 IE-02 
93  2  0.5000  -0.0500  0.1750 
85665.  -487,81  481.11  24702.  -82.856  -8328.7 

0.3095 1E-01-0.20127E-G1-0.26677E-02-0.44800E-02-0.98246E-04-0.98758E-02 


ELMDAT(  2)  =  2  for  Element  87 
CALCULATING  STRESSES 


ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 
87  1  0.5000  -0.1500  0.1250 

89013.  -2093.8  -402.32  -32834.  116.77  -559.10 

0.30985E-01-0.23030E-01-0.22899E-02  0.54083E-03  0.13846E-03-0.66295E-03 
87  2  0.5000  -0.1500  0.1750 

91559.  687.80  -355.99  35803.  -49.614  -563.02 

0.30953E-01-0.22922E-01-0.22899E-02-0.48002E-04-0.58830E-04-Q.66760E-03 


ELMDATf  2)  =  2  for  Element  81 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 
81  1  0.5000  -0.2500  0.1250 

92541.  511.04  46.429  -36300.  38.825  -45.634 

0.3 1 266E-0 1  -0.23296E-0 1 -0.20985E-02-0.96467E-04  0.45037E-04-0.54 1 1  IE-04 

81  2  0.5000  -0.2500  0.1750 

91228.  -503.58  43.513  34929.  43.803  -50.369 

0.3 1 176E-0 1-0.232 10E-0 1  -0.20985E-02-0. 1 5644E-03  0.5 1939E-04-0.59726E-04 


ELMDATf  2)  =  2  for  Element  75 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 
75  1  0.5000  -0.3500  0.1250 

93056.  1234.2  69.668  -36795.  -24.362  10.832 

0.31254E-01-Q.23186E-01-0.21 144E-02-0.10709E-03-0.28837E-04  0.12844E-04 

75  2  0.5000  -0.3500  0.1750 

92021.  192.20  60.987  35598.  -15.948  0.8C817 

0.3 1 249E-0 1 -0.23 195E-0 1  -0.2 1 144E-02-0, 1054  IE-03-0. 1891  IE-04  0.95829E-06 


ELMDAT(  2)  =  2  for  Element  69 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 
69  1  0.5000  -0.4500  0.1250 

92285.  732.56  48.222  -36100.  -12.680  12.178 

0.31 162E-0 1-0.23 1 17E-01-0.21 176E-02  0.401 14E-05-0.15035E-04  0.14440E-04 
69  2  0.5000  -0.4500  0.1750 

92654.  999.15  69.529  36388.  -18.344  1.2282 

0.3 1 2 10E-0 1  -0.23 1 30E-0 1  -0.2 1 1 76E-02  0. 20634E-04-0.2 175  IE-04  0. 145 64E-05 


ELMDAT(  2)  =  2  for  Element  63 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 
63  1  0.5000  -0.5500  0.1250 

91328.  -39.729  -12.260  -35374.  -29.173  9.7503 

0.3 1048E-0 1  -0.23 121 E-0 1  -0.2 1 125E-02  0.40684E-04-0.34592E-04  0. 1 1561E-04 
63  2  0.5000  -0.5500  0.1750 

92654.  1148.9  82.563  36468.  -27.817  2.3970 

0.31 167E-01-0.23085E-0 1-0.21 125E-02  0.33844E-O4-O.32984E-O4  0.28422E-05 


ELMDAT(  2)  =  2  for  Element  57 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 
57  1  0.5000  -0.6500  0.1250 

90467.  -266.98  -471.66  -34788.  -3.1402  438.84 

C.30879E-0 1  -0.229 16E-0 1-0.23220E-02  0.18055E-03-0.37235E-05  0.52036E-03 

57  2  0.5000  -0.6500  0.1750 

91468.  411.51  -311.82  35112.  172.14  443.69 

0.3 1 1G5E-0 1-0.22880E-0 1 -0.23220E-02-0. 340 1 0E-O3  0.20412E-03  0.52610E-03 


ELMDAT(  ?.)=  2  for  Element  51 
CALCULATING  STRESSES 

ELEM  X  Y  Z 

SX  SY  SZ  SXY  SYZ  SXZ 
51  1  0.5000  -0.7500  0.1250 

85291.  -459.24  51.317  -25167.  -33.788  8122.8 

0.307 19E-0 1  -0.20 1 20E-01  -0.27950E-02  0.42034E-02-0.40064E-04  0.96317E-02 
51  2  0.5000  -0.7500  0.1750 
86481.  279.66  250.74  25474.  189.82  8152.0 

0.3 1 0 1 6E-0 1  -0.2009 1 E-0 1 -0. 27950E-02-0.44209E-02  0. 22508E-03  0.96663E-02 


